## # A tibble: 19 x 5 ## dist successes tries p error_sd ## <int> <dbl> <dbl> <dbl> <dbl> ## 1 2 1346 1443 0.9327789 0.006591874 ## 2 3 577 694 0.8314121 0.014211556 ## 3 4 337 455 0.7406593 0.020546548 ## 4 5 208 353 0.5892351 0.026185062 ## 5 6 149 272 0.5477941 0.030178131 ## 6 7 136 256 0.5312500 0.031188905 ## 7 8 111 240 0.4625000 0.032183960 ## 8 9 69 217 0.3179724 0.031613007 ## 9 10 67 200 0.3350000 0.033374766 ## 10 11 75 237 0.3164557 0.030211036 ## 11 12 52 202 0.2574257 0.030762402 ## 12 13 46 192 0.2395833 0.030803744 ## 13 14 54 174 0.3103448 0.035072250 ## 14 15 28 167 0.1676647 0.028907578 ## 15 16 27 201 0.1343284 0.024052622 ## 16 17 31 195 0.1589744 0.026184896 ## 17 18 33 191 0.1727749 0.027354921 ## 18 19 20 147 0.1360544 0.028277490 ## 19 20 24 152 0.1578947 0.029576394
06/02/2018
Example Data - Golfing
Visualise the data
Logistic Regression?
logistic_binom <- glm(p ~ dist,
weights = tries,
data = golf,
family = binomial(link = "logit"))
logistic_binom
## ## Call: glm(formula = p ~ dist, family = binomial(link = "logit"), data = golf, ## weights = tries) ## ## Coefficients: ## (Intercept) dist ## 2.2312 -0.2557 ## ## Degrees of Freedom: 18 Total (i.e. Null); 17 Residual ## Null Deviance: 2411 ## Residual Deviance: 255.3 AIC: 365.9
… in Stan (the easy version)
library(rstanarm)
stan_logistic_binom <- stan_glm(p ~ dist,
weights = tries,
data = golf,
family = binomial(link = "logit"))
## ## SAMPLING FOR MODEL 'binomial' NOW (CHAIN 1). ## ## Gradient evaluation took 9.8e-05 seconds ## 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds. ## Adjust your expectations accordingly! ## ## ## Iteration: 1 / 2000 [ 0%] (Warmup) ## Iteration: 200 / 2000 [ 10%] (Warmup) ## Iteration: 400 / 2000 [ 20%] (Warmup) ## Iteration: 600 / 2000 [ 30%] (Warmup) ## Iteration: 800 / 2000 [ 40%] (Warmup) ## Iteration: 1000 / 2000 [ 50%] (Warmup) ## Iteration: 1001 / 2000 [ 50%] (Sampling) ## Iteration: 1200 / 2000 [ 60%] (Sampling) ## Iteration: 1400 / 2000 [ 70%] (Sampling) ## Iteration: 1600 / 2000 [ 80%] (Sampling) ## Iteration: 1800 / 2000 [ 90%] (Sampling) ## Iteration: 2000 / 2000 [100%] (Sampling) ## ## Elapsed Time: 0.110557 seconds (Warm-up) ## 0.114992 seconds (Sampling) ## 0.225549 seconds (Total) ## ## ## SAMPLING FOR MODEL 'binomial' NOW (CHAIN 2). ## ## Gradient evaluation took 2e-05 seconds ## 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds. ## Adjust your expectations accordingly! ## ## ## Iteration: 1 / 2000 [ 0%] (Warmup) ## Iteration: 200 / 2000 [ 10%] (Warmup) ## Iteration: 400 / 2000 [ 20%] (Warmup) ## Iteration: 600 / 2000 [ 30%] (Warmup) ## Iteration: 800 / 2000 [ 40%] (Warmup) ## Iteration: 1000 / 2000 [ 50%] (Warmup) ## Iteration: 1001 / 2000 [ 50%] (Sampling) ## Iteration: 1200 / 2000 [ 60%] (Sampling) ## Iteration: 1400 / 2000 [ 70%] (Sampling) ## Iteration: 1600 / 2000 [ 80%] (Sampling) ## Iteration: 1800 / 2000 [ 90%] (Sampling) ## Iteration: 2000 / 2000 [100%] (Sampling) ## ## Elapsed Time: 0.112554 seconds (Warm-up) ## 0.105586 seconds (Sampling) ## 0.21814 seconds (Total) ## ## ## SAMPLING FOR MODEL 'binomial' NOW (CHAIN 3). ## ## Gradient evaluation took 2.2e-05 seconds ## 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds. ## Adjust your expectations accordingly! ## ## ## Iteration: 1 / 2000 [ 0%] (Warmup) ## Iteration: 200 / 2000 [ 10%] (Warmup) ## Iteration: 400 / 2000 [ 20%] (Warmup) ## Iteration: 600 / 2000 [ 30%] (Warmup) ## Iteration: 800 / 2000 [ 40%] (Warmup) ## Iteration: 1000 / 2000 [ 50%] (Warmup) ## Iteration: 1001 / 2000 [ 50%] (Sampling) ## Iteration: 1200 / 2000 [ 60%] (Sampling) ## Iteration: 1400 / 2000 [ 70%] (Sampling) ## Iteration: 1600 / 2000 [ 80%] (Sampling) ## Iteration: 1800 / 2000 [ 90%] (Sampling) ## Iteration: 2000 / 2000 [100%] (Sampling) ## ## Elapsed Time: 0.100435 seconds (Warm-up) ## 0.107553 seconds (Sampling) ## 0.207988 seconds (Total) ## ## ## SAMPLING FOR MODEL 'binomial' NOW (CHAIN 4). ## ## Gradient evaluation took 2.2e-05 seconds ## 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds. ## Adjust your expectations accordingly! ## ## ## Iteration: 1 / 2000 [ 0%] (Warmup) ## Iteration: 200 / 2000 [ 10%] (Warmup) ## Iteration: 400 / 2000 [ 20%] (Warmup) ## Iteration: 600 / 2000 [ 30%] (Warmup) ## Iteration: 800 / 2000 [ 40%] (Warmup) ## Iteration: 1000 / 2000 [ 50%] (Warmup) ## Iteration: 1001 / 2000 [ 50%] (Sampling) ## Iteration: 1200 / 2000 [ 60%] (Sampling) ## Iteration: 1400 / 2000 [ 70%] (Sampling) ## Iteration: 1600 / 2000 [ 80%] (Sampling) ## Iteration: 1800 / 2000 [ 90%] (Sampling) ## Iteration: 2000 / 2000 [100%] (Sampling) ## ## Elapsed Time: 0.113862 seconds (Warm-up) ## 0.090601 seconds (Sampling) ## 0.204463 seconds (Total)
… in Stan
coefficients(logistic_binom)
## (Intercept) dist ## 2.2312106 -0.2556919
coefficients(stan_logistic_binom)
## (Intercept) dist ## 2.2342119 -0.2560822
Priors
prior_summary(stan_logistic_binom)
## Priors for model 'stan_logistic_binom'
## ------
## Intercept (after predictors centered)
## ~ normal(location = 0, scale = 10)
##
## Coefficients
## ~ normal(location = 0, scale = 2.5)
## **adjusted scale = 0.44
## ------
## See help('prior_summary.stanreg') for more details
stan_logistic_binom2 <- stan_glm(p ~ dist,
weights = tries,
data = golf,
family = binomial(link = "logit"),
prior = laplace(location = 0, scale = 10),
prior_intercept = cauchy(location = 0, scale = 2.5))
Predictions
predicted_proportions <- predict(logistic_binom, type='response') str(predicted_proportions)
## Named num [1:19] 0.848 0.812 0.77 0.722 0.668 ... ## - attr(*, "names")= chr [1:19] "1" "2" "3" "4" ...
stan_linepred <- posterior_linpred(stan_logistic_binom, transform = T) str(stan_linepred)
## num [1:4000, 1:19] 0.849 0.846 0.839 0.848 0.842 ... ## - attr(*, "dimnames")=List of 2 ## ..$ iterations: NULL ## ..$ : chr [1:19] "1" "2" "3" "4" ...
Posterior Predictive
logist_binom_predict <- posterior_predict(stan_logistic_binom, draws = 500) str(logist_binom_predict)
## 'ppd' int [1:500, 1:19] 1231 1220 1196 1224 1245 1194 1261 1226 1224 1226 ... ## - attr(*, "dimnames")=List of 2 ## ..$ : NULL ## ..$ : chr [1:19] "1" "2" "3" "4" ...
# convert sucesses to proportions logist_binom_predict <- t(t(logist_binom_predict) / golf$tries) str(logist_binom_predict)
## ppd [1:500, 1:19] 0.853 0.845 0.829 0.848 0.863 ... ## - attr(*, "dimnames")=List of 2 ## ..$ : NULL ## ..$ : chr [1:19] "1" "2" "3" "4" ...
\[ p(y^{rep}|y) = \int_{\Theta} p(y^{rep}|\theta) p(\theta|y)d\theta \]
A Geometric Model
\[ P(\text{success from dist } x) \\ = 2\Phi\Bigg(\frac{1}{\sigma} arcsin\Big(\frac{R-r}{x}\Big)\Bigg)-1 \]
P_success <- function(sigma, R, r, x){
2 * pnorm(asin((R-r)/x) / sigma) - 1
}
x <- seq(0,30,0.1)
qplot(x, P_success(1*pi/180, R, r, x),
geom="line") + ylim(0,1)
Stan Model Part 1
data {
int N;
int<lower=0> tries[N];
int<lower=0> successes[N];
real<lower=0> dist[N];
real R;
real r;
}
parameters {
real<lower = 0> sigma;
}
Stan Model Part 2
model {
real p[N];
for (n in 1:N)
p[n] = 2 * Phi(asin((R - r) / dist[n]) / sigma) - 1;
successes ~ binomial(tries, p);
}
Vectorised:
model {
real p[N];
p = 2 * Phi(asin((R - r) / dist) / sigma) - 1;
// A comment
successes ~ binomial(tries, p);
}
Stan Model Part 3
generated quantities {
vector[N] sucess_predictions;
for (n in 1:N) {
real p;
p = 2 * Phi(asin((R - r) / dist[n]) / sigma) - 1;
sucess_predictions[n] = binomial_rng(tries[n], p);
}
}
Stan Data Types
int N; real x; vector[N] x; simplex[K] x; unit_vector[K] x; ordered[K] x; positive_ordered[K] x; cov_matrix[K] x; corr_matrix[K] x; cholesky_factor_cov[K] x; cholesky_factor_corr[K] x;
Prepare Data & Compile Model
str(data)
## List of 6 ## $ N : num 19 ## $ tries : num [1:19] 1443 694 455 353 272 ... ## $ successes: num [1:19] 1346 577 337 208 149 ... ## $ dist : int [1:19] 2 3 4 5 6 7 8 9 10 11 ... ## $ r : num 0.07 ## $ R : num 0.177
library(rstan) #options(mc.cores = parallel::detectCores()) golf_stan <- stan_model(file = "golf.stan")
## In file included from filecd95fc1d82f.cpp:8: ## In file included from /Users/wells/R_libraries/StanHeaders/include/src/stan/model/model_header.hpp:4: ## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math.hpp:4: ## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math/rev/mat.hpp:4: ## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math/rev/core.hpp:12: ## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5: ## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math/rev/core/var.hpp:7: ## In file included from /Users/wells/R_libraries/BH/include/boost/math/tools/config.hpp:13: ## In file included from /Users/wells/R_libraries/BH/include/boost/config.hpp:39: ## /Users/wells/R_libraries/BH/include/boost/config/compiler/clang.hpp:200:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined] ## # define BOOST_NO_CXX11_RVALUE_REFERENCES ## ^ ## <command line>:6:9: note: previous definition is here ## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1 ## ^ ## 1 warning generated.
Sample!
golf_fit <- sampling(golf_stan, data=data)
## ## SAMPLING FOR MODEL 'golf' NOW (CHAIN 1). ## ## Gradient evaluation took 2e-05 seconds ## 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds. ## Adjust your expectations accordingly! ## ## ## Iteration: 1 / 2000 [ 0%] (Warmup) ## Iteration: 200 / 2000 [ 10%] (Warmup) ## Iteration: 400 / 2000 [ 20%] (Warmup) ## Iteration: 600 / 2000 [ 30%] (Warmup) ## Iteration: 800 / 2000 [ 40%] (Warmup) ## Iteration: 1000 / 2000 [ 50%] (Warmup) ## Iteration: 1001 / 2000 [ 50%] (Sampling) ## Iteration: 1200 / 2000 [ 60%] (Sampling) ## Iteration: 1400 / 2000 [ 70%] (Sampling) ## Iteration: 1600 / 2000 [ 80%] (Sampling) ## Iteration: 1800 / 2000 [ 90%] (Sampling) ## Iteration: 2000 / 2000 [100%] (Sampling) ## ## Elapsed Time: 0.025384 seconds (Warm-up) ## 0.025904 seconds (Sampling) ## 0.051288 seconds (Total) ## ## ## SAMPLING FOR MODEL 'golf' NOW (CHAIN 2). ## ## Gradient evaluation took 7e-06 seconds ## 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. ## Adjust your expectations accordingly! ## ## ## Iteration: 1 / 2000 [ 0%] (Warmup) ## Iteration: 200 / 2000 [ 10%] (Warmup) ## Iteration: 400 / 2000 [ 20%] (Warmup) ## Iteration: 600 / 2000 [ 30%] (Warmup) ## Iteration: 800 / 2000 [ 40%] (Warmup) ## Iteration: 1000 / 2000 [ 50%] (Warmup) ## Iteration: 1001 / 2000 [ 50%] (Sampling) ## Iteration: 1200 / 2000 [ 60%] (Sampling) ## Iteration: 1400 / 2000 [ 70%] (Sampling) ## Iteration: 1600 / 2000 [ 80%] (Sampling) ## Iteration: 1800 / 2000 [ 90%] (Sampling) ## Iteration: 2000 / 2000 [100%] (Sampling) ## ## Elapsed Time: 0.025744 seconds (Warm-up) ## 0.026367 seconds (Sampling) ## 0.052111 seconds (Total) ## ## ## SAMPLING FOR MODEL 'golf' NOW (CHAIN 3). ## ## Gradient evaluation took 8e-06 seconds ## 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds. ## Adjust your expectations accordingly! ## ## ## Iteration: 1 / 2000 [ 0%] (Warmup) ## Iteration: 200 / 2000 [ 10%] (Warmup) ## Iteration: 400 / 2000 [ 20%] (Warmup) ## Iteration: 600 / 2000 [ 30%] (Warmup) ## Iteration: 800 / 2000 [ 40%] (Warmup) ## Iteration: 1000 / 2000 [ 50%] (Warmup) ## Iteration: 1001 / 2000 [ 50%] (Sampling) ## Iteration: 1200 / 2000 [ 60%] (Sampling) ## Iteration: 1400 / 2000 [ 70%] (Sampling) ## Iteration: 1600 / 2000 [ 80%] (Sampling) ## Iteration: 1800 / 2000 [ 90%] (Sampling) ## Iteration: 2000 / 2000 [100%] (Sampling) ## ## Elapsed Time: 0.025561 seconds (Warm-up) ## 0.024421 seconds (Sampling) ## 0.049982 seconds (Total) ## ## ## SAMPLING FOR MODEL 'golf' NOW (CHAIN 4). ## ## Gradient evaluation took 7e-06 seconds ## 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. ## Adjust your expectations accordingly! ## ## ## Iteration: 1 / 2000 [ 0%] (Warmup) ## Iteration: 200 / 2000 [ 10%] (Warmup) ## Iteration: 400 / 2000 [ 20%] (Warmup) ## Iteration: 600 / 2000 [ 30%] (Warmup) ## Iteration: 800 / 2000 [ 40%] (Warmup) ## Iteration: 1000 / 2000 [ 50%] (Warmup) ## Iteration: 1001 / 2000 [ 50%] (Sampling) ## Iteration: 1200 / 2000 [ 60%] (Sampling) ## Iteration: 1400 / 2000 [ 70%] (Sampling) ## Iteration: 1600 / 2000 [ 80%] (Sampling) ## Iteration: 1800 / 2000 [ 90%] (Sampling) ## Iteration: 2000 / 2000 [100%] (Sampling) ## ## Elapsed Time: 0.025476 seconds (Warm-up) ## 0.024774 seconds (Sampling) ## 0.05025 seconds (Total)
Sampling Options
golf_fit <- stan(file = "golf.stan",
model_name = "golf_putting",
data = data,
chains = 3,
iter = 1000,
warmup = 200,
thin = 2,
seed = 42,
cores = 2,
algorithm = "HMC")
golf_fit <- vb(golf_stan, data=data)
make golf ./golf sample data file=golf_data.R ./golf variational data file=golf_data.R
print(golf_fit)
## Inference for Stan model: golf. ## 4 chains, each with iter=2000; warmup=1000; thin=1; ## post-warmup draws per chain=1000, total post-warmup draws=4000. ## ## mean se_mean sd 2.5% 25% 50% ## sigma 0.03 0.00 0.00 0.03 0.03 0.03 ## sucess_predictions[1] 1378.59 0.17 9.00 1360.00 1373.00 1379.00 ## sucess_predictions[2] 568.41 0.19 11.08 547.00 561.00 568.00 ## sucess_predictions[3] 311.62 0.18 10.37 291.00 305.00 312.00 ## sucess_predictions[4] 203.93 0.16 9.76 185.00 197.00 204.00 ## sucess_predictions[5] 135.01 0.14 8.34 119.00 129.00 135.00 ## sucess_predictions[6] 111.08 0.13 8.16 95.00 106.00 111.00 ## sucess_predictions[7] 92.38 0.13 7.63 77.00 87.00 92.00 ## sucess_predictions[8] 74.92 0.11 6.98 62.00 70.00 75.00 ## sucess_predictions[9] 62.37 0.11 6.70 50.00 58.00 62.00 ## sucess_predictions[10] 67.47 0.12 6.98 54.00 63.00 67.00 ## sucess_predictions[11] 52.98 0.10 6.31 41.00 49.00 53.00 ## sucess_predictions[12] 46.54 0.10 5.99 35.00 43.00 47.00 ## sucess_predictions[13] 39.26 0.09 5.50 29.00 36.00 39.00 ## sucess_predictions[14] 35.16 0.08 5.28 25.00 32.00 35.00 ## sucess_predictions[15] 39.88 0.09 5.61 29.00 36.00 40.00 ## sucess_predictions[16] 36.41 0.09 5.40 26.00 33.00 36.00 ## sucess_predictions[17] 33.74 0.08 5.28 24.00 30.00 34.00 ## sucess_predictions[18] 24.74 0.07 4.53 16.00 22.00 25.00 ## sucess_predictions[19] 24.14 0.07 4.50 16.00 21.00 24.00 ## lp__ -2926.77 0.02 0.72 -2928.80 -2926.92 -2926.48 ## 75% 97.5% n_eff Rhat ## sigma 0.03 0.03 1719 1 ## sucess_predictions[1] 1385.00 1396.00 2699 1 ## sucess_predictions[2] 576.00 590.00 3270 1 ## sucess_predictions[3] 319.00 331.03 3227 1 ## sucess_predictions[4] 210.25 223.00 3691 1 ## sucess_predictions[5] 141.00 151.00 3473 1 ## sucess_predictions[6] 116.00 127.00 3790 1 ## sucess_predictions[7] 97.00 107.00 3652 1 ## sucess_predictions[8] 80.00 89.00 4000 1 ## sucess_predictions[9] 67.00 76.00 4000 1 ## sucess_predictions[10] 72.00 81.00 3678 1 ## sucess_predictions[11] 57.00 65.00 3957 1 ## sucess_predictions[12] 50.00 59.00 3953 1 ## sucess_predictions[13] 43.00 50.00 4000 1 ## sucess_predictions[14] 39.00 46.00 4000 1 ## sucess_predictions[15] 44.00 51.00 4000 1 ## sucess_predictions[16] 40.00 47.00 3681 1 ## sucess_predictions[17] 37.00 44.00 4000 1 ## sucess_predictions[18] 28.00 34.00 4000 1 ## sucess_predictions[19] 27.00 33.00 3944 1 ## lp__ -2926.31 -2926.27 1671 1 ## ## Samples were drawn using NUTS(diag_e) at Sun Feb 4 00:56:23 2018. ## For each parameter, n_eff is a crude measure of effective sample size, ## and Rhat is the potential scale reduction factor on split chains (at ## convergence, Rhat=1).
Diagnostics
stan_trace(golf_fit, pars = "sigma")
stan_ac(golf_fit, pars = 'sigma')
Extract Samples
str(as.matrix(golf_fit))
## num [1:4000, 1:21] 0.0267 0.0266 0.0265 0.0266 0.0267 ... ## - attr(*, "dimnames")=List of 2 ## ..$ iterations: NULL ## ..$ parameters: chr [1:21] "sigma" "sucess_predictions[1]" "sucess_predictions[2]" "sucess_predictions[3]" ...
str(as.array(golf_fit))
## num [1:1000, 1:4, 1:21] 0.0267 0.0266 0.0265 0.0266 0.0267 ... ## - attr(*, "dimnames")=List of 3 ## ..$ iterations: NULL ## ..$ chains : chr [1:4] "chain:1" "chain:2" "chain:3" "chain:4" ## ..$ parameters: chr [1:21] "sigma" "sucess_predictions[1]" "sucess_predictions[2]" "sucess_predictions[3]" ...
str(rstan::extract(golf_fit))
## List of 3 ## $ sigma : num [1:4000(1d)] 0.027 0.027 0.0261 0.0272 0.0273 ... ## ..- attr(*, "dimnames")=List of 1 ## .. ..$ iterations: NULL ## $ sucess_predictions: num [1:4000, 1:19] 1364 1382 1384 1372 1363 ... ## ..- attr(*, "dimnames")=List of 2 ## .. ..$ iterations: NULL ## .. ..$ : NULL ## $ lp__ : num [1:4000(1d)] -2927 -2927 -2927 -2927 -2927 ... ## ..- attr(*, "dimnames")=List of 1 ## .. ..$ iterations: NULL
Check Fit
Posterior Predictive
One liner PPCs with library(bayesplot)
ppc_violin_grouped(y=golf$p, yrep=logist_binom_predict[1:500,], group = golf$dist, y_size = 2, y_draw = "both") + theme(legend.position = "bottom")
ppc_violin_grouped(y=golf$p, yrep=golf_predictive[1:500,], group = golf$dist, y_size = 2, y_draw = "both") + theme(legend.position = "bottom")
ppc_stat_grouped(golf$p, logist_binom_predict[1:500,], group = golf$dist, stat = "median") + theme(legend.position = "bottom")
ppc_stat_grouped(golf$p, golf_predictive[1:500,], group = golf$dist, stat = "median") + theme(legend.position = "bottom")
mcmc_areas(as.matrix(golf_fit)[,c(11,12,13), drop=FALSE])
ppc_dens_overlay(golf$p,
golf_predictive[1:500,])
ppc_dens_overlay(golf$p,
logist_binom_predict[1:500,])
Shinystan
library(shinystan) launch_shinystan(golf_fit)
Alternative 'Sampling' Statements
beta ~ normal(0, 1); target += normal_lpdf(beta | 0, 1); target += -0.5 * beta^2; target += -0.5 * beta’ * beta; target += -0.5 * dot_self(beta);
Stan Code Linear Regression
data {
int N ; # integer, number of observations
int K ; # integer, number of columns in model matrix
matrix[N,K] X ; # N by K model matrix
vector[N] y ; # vector of N observations
}
parameters {
real<lower=0> sigma ; # real number > 0, standard deviation
vector[K] beta ; # K-vector of regression coefficients
}
model {
beta ~ normal(0, 5) ; # prior for betas
sigma ~ cauchy(0, 2.5) ; # prior for sigma
y ~ normal(X*beta, sigma) ; # vectorized likelihood
}
generated quantities {
vector[N] y_rep ; # vector of same length as the data y
for (n in 1:N)
y_rep[n] <- normal_rng(X[n]*beta, sigma) ;
}